Hole Localization in Molecular Crystals From Hybrid Density Functional Theory 



and Paul F. Barbargfl] 

Center for Nano and Molecular Science and Technology, 
The University of Texas at Austin, Austin, Texas, 78712, USA 

Kevin Leun^ 

Surface and Interface Sciences Department, MS 14 15, 
Sandia National Laboratory, Albuquerque, New Mexico 87185, USA 
(Dated; January 20, 2013) 

We use first-principles computational methods to examine hole trapping in organic molecular 
crystals. We present a computational scheme based on the tuning of the fraction of exact exchange 
in hybrid density functional theory to eliminate the many-electron self-interaction error. With small 
organic molecules, we show that this scheme gives accurate descriptions of ionization and dimer 
dissociation. We demonstrate that the excess hole in perfect molecular crystals form self-trapped 
molecular polarons. The predicted absolute ionization potentials of both localized and delocalized 
holes are consistent with experimental values. 



Charge localization in organic solids and interfaces has 
important implications for the functioning of organic de- 
vices [U 12 . It has been long suggested that charge carri- 
ers in organic molecular crystals may self-trap to form lo- 
calized small polarons [1-4^ through interaction with the 
surrounding electron and lattice. Models based on the 
concept of polaron have led to theoretical understanding 
of charge transport of these materials (see e.g. [HH]). On 
the other hand, direct evidence and atomic scale probe of 
localized charge states in organic crystals remain lacking. 
Ab initio computations have recently been summoned to 
provide accurate quantitative description of polaron state 
in covalent solids [HI IZ] • In this work, we apply hybrid 
density functional theory (DFT) with carefully calibrated 
admixtures of exact exchange to molecular crystals and 
conclusively demonstrate from first principles the exis- 
tence of self-trapped hole polaron in defect-free molecular 
crystals. Our choice to work with small organic molecules 
permits systematic analysis of DFT exchange correlation 
(xc) functionals in both gas and solid phases. It also 
allows predictions of the absolute ionization potentials 
(IP) of both polaronic localized and free delocalized hole 
state in solids, revealing large errors in semilocal (si) DFT 
functionals that to our knowledge have never been dis- 
cussed in the literature. This successfully treatment of 
hole localization and delocalization properties will aid the 
understanding of charge trapping [51 H] and carrier mo- 
bility [lOj that are important for organic electronics. 

A key prerequisite to such studies is the correct choice 
of DFT functionals. Previous theoretical studies reveal 
an extreme sensitivity of charge localization to theoret- 
ical approximations [SI [3 [TTmi] . Standard semilocal 
functionals for the xc energy, such as the local density 
approximation (LDA) and generalized gradient approxi- 
mation (GGA), have been successful in predicting atomic 
structures and electronic properties of many molecular 
complexes and solid state materials. But they tend to 
fail in systems involving fractional charges, e.g. they dis- 



sociate molecular ion into 2 H°-^+ rather than into 
H''" and H. The problem stems from self-interaction error 
(SIE) in the standard density functionals [15 . For a one 
electron system SIE is conveniently described as the in- 
complete cancellation of spurious self-exchange and self- 
Coulomb energy [16] . A more general many-electron SIE 
(also known as the delocalization error) [T71 [T5] arises 
from the incorrect behavior of the total energy as a func- 
tional of the average number of electrons. While the ex- 
act energy varies linearly between adjacent integers jl9) . 
commonly used xc functionals give rise to convex or con- 
cave behavior [2Dj, which tends to overstabilize either 
the delocalized or the localized electron density respec- 
tively. This will clearly impair accurate prediction of 
charge trapping and localization in organic solids. To 
achieve a many-electron-SIE-free situation, the ground 
state total energy E{N) should be linear in the electron 
number N and has a realistic slope [H]. 

In this letter, we achieve this condition by tuning the 
mixing parameter a in a hybrid functional, which com- 
bines a fraction of the nonlocal exact Hartree-Fock (HE) 
type exchange with the semilocal exchange, 

= aE^^ + (1 - a)El' (1) 

where < a < 1 defines the mixing coefficient. The 
widely used PBEO hybrid functional [22l[23| sets a = 0.25 
and leads to band gaps and atomic energies in better 
agreement with the experiments. However, we will show 
that the global PBEO {a = 0.25), owing to the residual 
SIE, gives an erroneous dissociation limit in the benzene 
dimer cation and a delocalized charge in molecular solids 
with excess charge. By contrast, our hybrid functional 
leads to correct dissociative behavior of molecular dimer 
cations and yields single molecular polaron-like hole local- 
ization in molecular solids. We demonstrate the validity 
of this approach for a number of small organic molecules 
with a focus on polyacenes. 

We have performed the calculations in the VASP pro- 
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FIG. 1: Total energy of the benzene molecule as a function 
of the electron number A'^; the energy is zeroed out at A'^ = 0. 
a = and a — 0.25 correspond to the PBE and the global 
PBEO hybrid functional, and a = 1.0 corresponds to the full 
HF exchange with PBE correlations. Inset: the curvature 
of the total energy with respect to the electron number as 
a function of a for benzene (solid red), naphthalene (dashed 
green), and anthracene (doted blue). 



gram using projector-augmented wave potentials [Mll25j . 
All the calculations are spin unrestricted and without 
symmetry constraint. We use a kinetic cutoff of 400 
eV and T k-point Brillouin zone sampling. For charged 
systems we add the monopole correction [26]. Because 
the molecular solids considered in this work exhibit high 
symmetry and obvious physical boundaries (terminated 
as intact molecules), the large and unscreened quadruple 
correction ^27 can also be determined to yield absolute 
ionization potentials. 

First we consider isolated polyacene molecules. Fig. [l] 
illustrates the total energy E{N) of benzene as a func- 
tion of the electron number TV, where N varies from 
M - 1 (cation) to M (neutral). The PBE curve {a = 0) 
shows an incorrect convex behavior as the fraction of hole 
charge increases, and PBEO {a = 0.25) improves only 
slightly upon PBE but remains convex. As a increases, 
the energy curve becomes increasingly more linear while 
at a = 1.0, i.e., with the full HF exchange (plus PBE 
correlation), it becomes concave. The linear dependence 
of the curvature d'^E{N)/dN'^ on a (Fig.[l]inset) enables 
us to determine the critical value ac that fulfills the linear 
criterion. For all three polyacenes, the condition is satis- 
fied practically at the same value ac = 0.7 (much higher 
than a = 0.25 used in PBEO) suggesting no direct depen- 
dence of ac on the conjugation length. At a = 0.7, the 
gas phase ionization energies Ig = 9.13, 8.01, and 7.23 eV 
for benzene, naphthalene, and anthracene, in agreement 
with the respective experimental values of 9.24, 8.14, and 
7.45 eV dH]. The HOMO-LUMO gaps of 9.2, 8.01 and 
6.85 eV are also in agreement with the respective exper- 
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FIG. 2: Density of states of the benzene crystal 2x2x2 
supercell with an excess hole calculated from (a) PBE and (b) 
the hybrid functional at a = 0.7. The arrows guide the eyes 
for the hole states. The molecular distortion related states 
are visible in the gap in (b). 



imental values of 9.98, 8.0, and 6.76 eV pH] . 

For a benzene dimer cation (CeHg) J in the dissociation 
limit, we have found that both PBE and PBEO (even at 
a = 0.65) yield a symmetric dissociated state with -|-0.5e 
on each molecule. By contrast, the hybrid functional at 
a = 0.7 yields an asymmetric dissociated state, i.e., a 
neutral (CeHg)^ and a cation (CeHg)"*", as expected. 

With the protocol established in molecules, we proceed 
with the calculation of hole state in solids. Using a 2 x 2 x 
2 extended supercell consisting of 384 atoms and starting 
out with a perfect crystal at equilibrium, we insert an 
excess hole into the system by removing an electron and 
allow the atomic positions to fully relax while keeping the 
lattice parameters fixed at the experimental values [55] 

The density of states (DOS) for the positively charged 
benzene crystal is shown in Fig [2] The PBE hole state 
splits into half occupations at the top of the valence band 
leading to degenerate spin up and spin down states as 
shown in Fig. [2] (a). It is clear that this hole state is fully 
delocalized over the entire lattice. In the case of PBEO 
(a = 0.25), although the hole state lies well above the va- 
lence band, our calculation shows that it fully delocalizes 
in the crystal regardless of the structural relaxations. In 
contrast, the hybrid a = 0.7 functional, which exhibits 
no SIE for an isolated CeHg, predicts that the hole state 
occupies only a single spin state (see Fig. [2] (b)) which 
splits from the valence band by 3.0 eV. Upon relaxation, 
it yields a distorted lattice in which the spatial distribu- 
tion of the hole state is localized around a single molecule, 
as shown in Fig. [3] 

Next we analyze the ionization energies of the local- 
ized hole state {E^^^ — E^^^) and the delocalized hole 
i^dc\oc ~ -^?of)' where E+ and 5'°^^ refer to the energy 
of the positively charged and the neutral ground state 
equilibrium system. Table |l] shows the results for ben- 
zene and anthracene crystal. The delocalized hole energy 
is computed in the neutral equilibrium geometry. After 



TABLE I: Calculated ionization energy (in eV) of the local- 
ized and delocalized hole state and self-trapping energy of the 
polyacene crystals (E^^f is set zero). The absolute values of 
the ionization energy /loc and -^dcioc fl^re obtained by adding 
a quadruple correction (jjq which is independent of the extent 
of charge delocalization. Experimental value is taken from 
Ref. [33] 



polyacenes 




^Expt 


benzene 


2.53 3.2 5.32 7.85 8.52 -0.67 


7.58 


anthracene 


-0.33 0.74 6.13 5.80 6.87 -1.07 


5.70 



adding finite size corrections [31], we find Jioc ^ 7.85 eV 
for benzene and 5.8 eV for anthracene, in good agree- 
ment with the experimental 7.58 and 5.7 eV associated 
with the hole polaron state [U [33]. The localized hole 
states therefore exhibits a lower ionization potential by 
E^*- — /loc — /dcioc of —0.67 eV for benzene and —1.07 
eV for anthracene. This "self-trapping energy" exceed 
systematic uncertainties due to, say, electrostatic correc- 
tions [26] and clearly favors hole localization. We thus 
obtain a general prediction of intrinsic self-trapping of 
holes in a small polyacene crystals in the absence of im- 
purities and crystal defects. 

Compared to the gas phase, the stabilization energy for 
the delocalized hole in sohd A/dcioc = Ig — ^dcioc of 0.61 
eV and 0.36 eV for benzene and anthracene corresponds 
to the energy associated with spreading out the charge 
between adjacent molecules and is consistent with the 
calculated bandwidth ~ 0.5 eV of benzene and ~ 0.34 eV 
of anthracene crystal [34]. In contrast, the stabilization 
energy for the localized hole A/ioc of 1.28 eV and 1.43 eV 
for benzene and anthracene compares well with the total 
polarization energy measured experimentally [1] [53] and 
can be attributed to electronic polarization and geom- 
etry relaxation. The electronic polarization energy can 
be crudely estimated from a dielectric continuum model 
AP = -(1 - l/e)(e2/2i?) where R is the radius of the 
sphere surrounding the excess charge and e is the dielec- 
tric constant [T] . For benzene we estimate R on the order 
of 3-4 A (enclosing > 95% of the localized spin) from the 
spin distribution of the hole state in Fig. [3] which yields 
a polarization energy of ~ 1 — 1.4 eV in a range reported 
for polyacene crystals [TJ [35] . The geometry relaxation 
leads to an energy gain associated with the intramolecu- 
lar relaxation [37] of 0.17 eV for benzene and 0.14 eV for 
anthracene [35]. The latter value is slightly higher than 
previous obtained from B3LYP hybrid functional [36l [38] 
but compares well with that from the linear coupling con- 
stant method [35] . There is also a cost of energy due to 
elastic lattice distortion (i.e., E^^^ — E'^^^) of 0.36 eV for 
benzene and 0.15 eV for anthracene. The self-trapping we 
observed is thus a result of competition between energy 
gain and cost associated with localization and delocaliza- 
tion. The values of these energy components are similar 
to those identified for molecular polarons [Tl where the 
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FIG. 3: Spatial distribution of the excess hole state in the 
benzene crystal calculated from the hybrid a = 0.7 functional. 
The hole is shown self-trapped around a single benzene site. 



intramolecular relaxation energy corresponds to the for- 
mation energy and thus supporting the notion that the 
localized hole state is indeed a self-trapped hole polaron. 
Predicting the outcome is not straightforward, especially 
in the nonlinear near-critical regime since the different 
energy components are not necessarily additive. The del- 
icate balance highlights the importance of accurate DFT 
approaches in which the various energy contributions are 
properly accounted for. 

To investigate why the standard PBE and PBEO func- 
tional fail to predict hole localization, we have also calcu- 
lated the stabilization energy of the delocalized hole using 
the PBE and PBEO functional. We find A/doioc = 3.34 
and 2.37 for a = and 0.25, respectively, with a surpris- 
ingly strong dependence on a. (The IP values obtained 
for (C6H6)„ clusters, with n up to 12, are given in the 
Supporting Information and confirm the periodic bound- 
ary condition predictions.) Compared to 0.61 eV calcu- 
lated at a — 0.7, the values of A/dcioc associated with 
a — and 0.25 are much higher than the experimental 
total polarization energy in molecular crystals [U [33] and 
make the delocalized state unphysically favorable. This 
dovetails with our monomer-based analysis that a = 0.7 
is needed to eliminate the many-electron SIE — even 
though a completely delocalized hole exhibits a vanish- 
ing one-electron SIE. To our knowledge, such strong IP 
dependence on DFT functionals has not been reported 
before. The theoretical reasoning for the shift observed 
in the delocalized hole energy is related to the significant 
shift of the one-particle self energy for the valence band 
maximum in solids when DFT is replaced with the GW 
approximation |40j. 

We have examined other small organic molecules be- 
sides polyacenes, including p-benzoquinone (C6H4O2) 
and 2,2'-bithiophene (C8H6S2). The critical value of a 
is somewhat system dependent, e.g., it is 0.6 for bithio- 
phene and 0.5 for benzoquinone, all much greater than 
a = 0.25 used in the global PBEO functional. The de- 
pendence of ac on the type of molecule seems to sug- 
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gest specific differences in tlie response of local environ- 
ment, including the chemical species, molecular geome- 
tries, and bondings, upon ionization. We have identified 
similar self-trapped hole polaron on a single molecular 
site in the associated crystals indicating a general behav- 
ior of excess charge in small molecular solids. Finally we 
note that transport mechanisms in molecular crystals, in 
particular, high mobility organic crystals, remain under 
debate [H]. The ability of our approach to accurately 
describe the presence of polaron state in molecular crys- 
tals will be crucial for future studies of polaron dynamics 
purely from first principles. 

In summary, after showing that DFT delocalization 
error is detrimental for the study of excess charges in 
organic molecular system, we suggest that this problem 
can be overcome with a hybrid density functional the- 
ory that incorporates an appropriate fraction of exact 
exchange. Our scheme enables an accurate description 
of the ionization and dimer dissociation limit of small or- 
ganic molecules and predicts self-trapped hole polaron in 
a perfect molecular crystal that is not captured in the 
standard DFT. Our result suggests that semi-local and 
even standard hybrid DFT functionals overestimate the 
tendency of hole delocalization owing to errors in the 
delocalized state. The gas phase ionization potentials of 
small molecules, arguably related to the energies of local- 
ized holes in solids, are however accurately predicted by 
even the semi-local functionals that we have considered. 
Our study may contribute to the understanding of pos- 
sible mechanisms of deep-hole trap states |9j commonly 
observed in organic electronics, and may potentially be 
applicable to covalently bonded solids [3 fTTHT^ . 
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